%--- main parameters
% SI System
%clear all;
function [] = Benchmark(option)
global Ux Uy  h  g  p  q r tmax dx dy Nx Ny y tmax;


dx1=1/dx;
dx2 = 1/(dx*dx);

dy1= 1/dy;
dy2 = 1/(dy*dy);

L=zeros(2*Nx*Ny,2*Nx*Ny);
A = generate_hyperbolic( Nx, Ny, Ux*0.5*dx1, Uy*0.5*dy1); 
L(1:Nx*Ny,1:Nx*Ny)=A;
L(Nx*Ny+1:2*Nx*Ny,Nx*Ny+1:2*Nx*Ny)=A;
A = generate_elliptic(Nx,Ny,h*dx2,h*dy2);
L(1:Nx*Ny,Nx*Ny+1:2*Nx*Ny)=A;

clear A;
for k=1:Nx*Ny
    L(Nx*Ny+k,k)=g;
end

 F = generate_elliptic(Nx,Ny,h*h/3*dx2,h*h/3*dy2);
 C=zeros(2*Nx*Ny,Nx*Ny);
 C(1:Nx*Ny,1:Nx*Ny)=F;
 
 E = generate_elliptic(Nx,Ny,-2/15*h*h*h*dx2,-2/15*h*h*h*dy2);
 
 for k=1:Nx*Ny
     E(k,k)=E(k,k)+h/3;
 end
 
 
 L = sparse(L);
 C=sparse(C);
 E=sparse(E);
 
dt=0.1;
 I = eye(2*Nx*Ny,2*Nx*Ny);
 Z1 = zeros(2*Nx*Ny,Nx*Ny);
 Z2 = zeros(Nx*Ny,2*Nx*Ny);
 Z3 = zeros(Nx*Ny,Nx*Ny);
 M = [I Z1; Z2 Z3];
tspan = 0:dt:tmax;  % time vector, 
%options1=odeset('RelTol',1e-8,'Stats','on'); 
 options2=odeset('RelTol',1e-6,'Stats','on', 'Mass', M,'MassSingular', 'yes'); 

tic

y0=[p;q;r];

%y0 = [p;q];

if option ==1
    disp('ODE 45 Explicit Method Selected');
[t,y]=ode45(@waves,tspan,y0,options1);
else
    disp('ODE 15s Implicit Method Selected');
   % [t,y]=ode15s(@waves,tspan,y0,options2, L,C,E);
    [t,y]=ode15s(@(t,y) waves(t,y,L,C,E),tspan,y0,options2);

    %[t,y]=ode15s(@upwinding,tspan,y0,options2);

end
toc


